Electronic Stopping and Momentum Density of Diamond Obtained from 

First-Principles Calculations 
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We calculate the "head" element or the (0, 0)-element of the wave-vector and frequency-dependent 
dielectric matrix of bulk crystals via first-principles, all-electron Kohn-Sham states in the integral of 
the irreducible polarizability in the random phase approximation. We approximate the macroscopic 
"head" element of the inverse matrix by its reciprocal value, and integrate over frequencies and 
momenta to obtain the electronic energy loss of protons at low velocities. Numerical evaluation for 
diamond targets predicts that the band gap causes a strong non-linear reduction of the electronic 
stopping power at ion velocities below 0.2 atomic units. 



PACS numbers: 61.80.Jh, 34.50.Bw, 71.45.Gm 
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I. FORMULATION AND METHOD 

The objective of this research is a quantitative, first- 
principles description of the energy deposition by a bare 
ion in diamond. The energy loss per unit path length 
of a massive, punctiformj charged particle with charge 
number Z\ to a target istH 
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dx 



v) 



d 3 q 



v q 



du>S(u> — q ■ v)lmK (q, u> 
(1) 



with v the projectile velocity in the target rest frame, 
e the elementary charge unit, and eo the vacuum per- 
mittivity. For reciprocal lattice vectors G, K is the 
G = G' = ("head") component of the inverse dielec- 
tric matrix defined by 



G' 



£ G ,G'(q, w)-KG',G"(q, W ) = ^GG" 



with respect to £G,G'(k, uj), the wave- vector and fre- 
quency dependent microscopic dielectric matrix. In 
terms of the irreducible polarizability II, the dielectric 
matrix is 



e G ,G'(q, u) = 6 G g> 



eo|q + G| 



r n(q+G,q+G',^), 



given in the Random Phase Approximationa (RPA) as a 
sum over all band pairs (z/, v 1 ) and an integral over the 
first Brillouin zone (BZ) 
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E v ^ are the band energies, and f u ^ are Fermi occupation 
numbers (/ y ,k = or 2 for E v ^ above or below the Fermi 
energy Ep). The matrix elements in Eq. (0) are 



(z/k 
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:(G+q)-r 



tdc). 



We determine the ground-state electronic structure of 
solids within Density Functional Theory (DFT) as es- 
tablished in the Kohn-Sham (KS) variational procedun 
and implemented in the computational package GTOFF 
The results of the all-electron, full-potential calculations 
are eigcnfunctions ip Ut k(r), expressed as linear combina- 
tion of Gaussian Type Orbitals (GTO's), and eigenvalues 

-Ei/,k- 

The Bloch functions are finally expanded in a (trun- 
'cated) plane wave (PW) series, 
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to represent tog as a simple sum over products of ex- 
pansion coefficients u^ cp where Vuc is the volume of 
a unit cell (UC). The equivalence of the GTO and PW 
representations is maintained by monitoring the accumu- 
lated norm for each \v, k) relative to the exact values, 



|M^k,c| 2 = Vac- 



(3) 
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We subdivide the integration region of the integral (|2j) 
into k-space tetrahedra. Recursive further subdivision 
of a given tetrahedron into smaller tetrahedra is done if 
/i/k — jvk+q is not constant over all four vertices. Next 
comes linearization of the product of the matrix elements 
in the numerator and of the energy denominator inside 
each tetrahedron for each to. Thejesulting approximated 
integral is evaluated analytically^ 

The product |q+ G| 2 eG,G'(q, u>) is calculated and, in 
compensation, the term q 2 in the denominator of Eq. 
(|l|) is dropped. The dielectric function is tabulated for 
q commensurate with the uniform mesh of wave vec- 
tors used in the underlying GTOFF calculation but cov- 
ering higher BZ's as well as the first. [Values outside 
the (q, w)-meshes are replaced by the vacuum response, 
IiTLKo.o(q, u>) = 0.] ImivG,G' is linearized inside each q- 
space tetrahedron. Multiplied by the linear factor q • v, 
the integrals (Q) over tetrahedra are done analytically, 
then summed. 
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II. DIAMOND 
A. Basis Sets, Electron Momentum Density 

Diamond is studied at the experimental lattice pa- 
rameter, a — 6.74071ao for the cubic unit cell, where 
ao denotes one _bohr. The Moruzzi-Janak- Williams 
parameterizationa of the Hedin-Lundqvist local-density 
approximation (LDA) to the exchange-correlation poten- 
tial is used. Small, highly contracted basis sets as used 
in Ref. |^ are generally insufficient to calculate the real 
parts of dielectric functions and subsequently the energw 
loss functions. We started from Partridge's 16sllp set,cl 
contracted the seven tightest s-functions and the four 
tightest p-functions, removed the most diffuse s and the 
two most diffuse p-functions to avoid approximate linear 
dependencies, and added a full set of three d-functions 
with exponents equal to the remaining most diffuse p- 
functions. (Without d-orbitals the total energy would 
rise by 0.24 eV/atom.) Site centered s- and /-type fit- 
ting functions were used with exponents for the s-types 
as in Ref. [?], and for the /-types 0.5/ and 0.2/ag as 
in Ref. ||. Space group and lattice type are those of sili- 
con; hence we may refer to a prior GTOFF studyla for the 
symmetry properties of the fitting functions. 

The density of states (DOS) computed with this 
9s6p3d basis (84 basis functions in the primitive UC) as 
shown in Fig. [j] is stable up to « 70 eV above the Fermi 
energy with respect to further de-contraction. Band gaps 
are too small compared with experiments, as usual for 
the LDAj-and known from other DFT calculations on 
diamond £3 

The lowest 24 bands were expanded into 531 plane 
waves, a cut-off at |G| = 7.3/ao, with a norm in Eq. (||) 
above 0.85Vuc for the A'-shell electrons (excitations from 
which were excluded in the subsequent calculations of 
£g.c arid dE/dx anyway), and a norm above 0.98Vuc for 
the remaining 22 bands. A first result is the all-electron 
momentum density (EMD) 

P (k + G)^-^^/,, k K, k , G | 2 

in Fig. ||. The values near zero momentum k + G = 
are the states with long wavelengths and reveal the 
macroscopic symmetry of the crystal system, the 4-fold 
axis of the cubic system here. If |k + G| is of the order 
of half a reciprocal lattice vector, the interference with 
the next nearest neighbors in the lattice becomes visible; 
the eight foothills in the figure may be interpreted as a 
projection of the four corners of a carbon tetrahedron 
onto the (001) plane complemented by the inversionE-H 

If |k + G| is large, the spherical symmetry of the core 
states prevails. 



B. Dielectric Function 

The element eo,o(q 7 ^) of the dielectric matrix is shown 
m Fig. |. The main absorption peak at 11 eV for |q| — > 
corresponds to direct transitions from the top of the 
valence to the bottom of the conduction band at X and 

lM 

Local Field Effects are illustrated in Fig. [|. Imeo.o 
repeats some values of Fig. ||; Imfl/A^o] includes an 
estimate of the local field by calculating Kg,g' as the 
inverse of a 9 x 9 dielectric matrix which contains G = 
and the eight vectors of the closest shell in the bcc recip- 
rocal lattice. The reduction of the values without LFE 
(|Imeo,o|, open symbols) compared to those with LFE 
(|Im(l / -ftTo.o) | , filled symbols) is of the order reported 
by Van Vechten and MartirO (without their "dynamical 
correlations"). The different sign of the effect for frequep. 
cies above and below the peak has been noticed before.Eil 
The differences are even smaller for the energy loss func- 
tion. Hence the energy loss in the next paragraph was 
calculated from eo,o(q, alone. 

C. Computed Electronic Energy Loss 

To integrate the energy loss function as described by 
Eq. (0), the BZ mesh was reduced to 8 x 8 x 8 points, 
i.e., 35 points in the irreducible Brillouin zone (IBZ). 
The result is shown in Fig. ^| The dielectric matrix 
was tabulated on a 30 x 30 x 30 mesh in q-space (par- 
allelepiped with three edges of length 6.1/ao), and the 
stopping power was integrated over the superset of all 
q-values obtained from these via point group operations 
with Seitz symbol {0|w}, 

*T G ,G'(q,w) = e lO ( G - G ') w if OG , OG ,(0q, W ). 

It is about 0-15% lower than experimental resultslll for 
v Rs 1 . . . 1.3wo which are shown in Fig. § {v « 2.19 • 10 6 
m/s is one atomic unit of velocity. Values within the 
LPDA represent intgrals of the all-electron density of 
GTOFF Fitting Functions weighted with the stopping 
number of the FEG in the RPA.) Reduction of the inte- 
grated q-region to the parallelepiped decreases the stop, 
ping cross section by approximately 10% at v « Vo.t3 
Within the framework of linear dielectric response an 
underestimation is appropriate though, because terms of 
O(Zf) will add about 15% to S(v) at v w 1.5w .O 

A new result is the "ionic" band gap, the non-linear 
suppression of S(v) for vS. 0.2vq, which is approximately 
the value extracted fronxJ 

«< Y\/ hLJ ° /E ° (4) 

for a band gap of huj g — 4 eV. (E is one rydberg.) 
By way of contrast, calculations within the local-plasma- 
density approximation usually integrate over volume el- 
ements in real space that are parameterized by the ho- 
mogeneous electron gas, and inevitably yield S(v) oc v at 



3 



low velocities. However, an attempt at experimental ver- 
ification of this reduction of the energy loss in a wide-gap 
material is handicapped by the additional nuclear energy 
loss, which has an estimated -maximum of S w 0.8 • 10~ 15 
eVcm 2 /atoms at v « 0.08vo-" 

Note that Eq. (Q) predicts contributions from K-sheW 
excitations to start at v » 2.2vq, which are not included 
here. They have been estimated to-shift the maximum 
to higher jitelocities by about 0.2t>oJi3 and to grow until 
v k 5.3wq-B3 
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FIG. 1: KS band-structure E v ^ using 16 x 16 x 16 k-points 
in the BZ (213 in the IBZ). The direct gap at T is 5.59 eV, the 
indirect gap 4.13 eV, and the width of the four valence bands 
with eight electrons 21.29 eV. The tjwp bands with four core 
electrons at —263 eV are not shownEj 




FIG. 2: The EMD as a function of k + G, which is varied 
in the plane spanned by (100) and (010). Both momentum 
components are measured in units of l/oo- Anisotropies of 
"projected" positron-EMD's are discussed in Ref. 
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FIG. 3: Calculated Ime G ,G' for G = G' = and q || (111). 
The underlying DFT calculation is based on a 12 x 12 x 12 
mesh of k-points in the BZ (98 points in the IBZ), which 
creates for this particular q-direction commensurate values 
at |q| = (j/12)27r\/3/a as shown here for j = 1, 2, 3, 
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FIG. 4: The absolute value of Imeo i o(q, ui) (open symbols) 
and of Iml/A'o,o(q, w) (filled symbols), q is parallel to (111) 
with qao = 0.135 (j = 1, squares) and gap = 0.269 (j = 2, 
circles) as in Fig. || Note that Refs. [u|[[4| and |2^ refer to the 
optical limit q = j = 0. 
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FIG. 5: Electronic stopping cross section of diamond. The 
stopping power was integrated on a mesh with 72 points on 
the tj-axis (0 . . . 103 eV). The lowest 28 bands were included 
in the sum overrband-pairs in Eq. (g). The line refers to the 
free electron gasEj with a density equivalent to 4 electrons per 
diamond atom. 
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FIG. 6: Experimental proton stopping r«ross section by 
Kaferbock et al.i3 and of graphite by Janni,l— I compared with 
calculated value*4or gaseous carbon and solid graphite from 
Kaneko's theory ,c3 and RPA free electron gas (FEG) valuesEJ 
with a homogeneous density equivalent to 4 or 6 electrons per 
diamond atom (Fermi velocity 1.457t>o or 1.668uo) as shown. 



